#--------------#
#### Prelim ####
#--------------#
library(tidyverse)
library(rdrobust)
library(stargazer)
library(foreign)
library(readstata13)
library(data.table)
library(xtable)
library(lfe)

# RDD output exporting functions:
source("code/rd.export.R")

#-------------#
#### Data  ####
#-------------#
sler4 <- read_csv("slers_processed.csv",guess_max = 10000)
ceda4 <- read_csv("ceda_processed.csv",guess_max = 61000)
mayors4 <- read_csv("mayors_processed.csv")

sler4 <- sler4 %>%
  mutate(lossmargin = 1-voteshare_toptwo-0.5)

ceda4 <- ceda4 %>%
  mutate(lossmargin = 1-voteshare_toptwo-0.5)

mayors4 <- mayors4 %>%
  mutate(lossmargin = 1-voteshare-0.5)

sler_firsttime <- sler4 %>% 
  group_by(candid) %>%
  arrange(year_month) %>%
  filter(row_number(year_month) == 1)

ceda4$elected <- ifelse(ceda4$voteshare_toptwo>=0.5,1,0)
ceda4$lost <- ifelse(ceda4$voteshare_toptwo<0.5,1,0)

ceda_firsttime <- ceda4 %>% 
  group_by(candid) %>%
  arrange(date) %>%
  filter(row_number(date) == 1)

mayors4 <- subset(mayors4, !duplicated(mayors4[,c("fips","YearData","month","name","lossmargin")]))

mayors_firsttime <- mayors4 %>% 
  group_by(candid) %>%
  arrange(YearData,month) %>%
  filter(row_number() == 1)

#----------------------#
#### SLERs analysis ####
#----------------------#
sler_all_cct <- with(subset(sler4), 
										 rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
sler_all_cct_90 <- with(subset(sler4), 
                     rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(sler_all_cct) # across all cands, effect of 0.312


sler_next4_all_cct <- with(subset(sler4), 
                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
summary(sler_next4_all_cct) # across all cands, effect of 0.341 on running in same office
sler_next4_all_cct_90 <- with(subset(sler4), 
                              rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))


sler_female_cct <- with(subset(sler4,female==1), 
												rdrobust(y=run_again_anywhere , x=lossmargin, c=0))
summary(sler_female_cct) # female cands, effect of 0.336
sler_female_cct_90 <- with(subset(sler4,female==1), 
                        rdrobust(y=run_again_anywhere , x=lossmargin, c=0, level=90))

sler_male_cct <- with(subset(sler4,female==0), 
												rdrobust(y=run_again_anywhere , x=lossmargin, c=0))
summary(sler_male_cct) # male cands, effect of 0.304
sler_male_cct_90 <- with(subset(sler4,female==0), 
                      rdrobust(y=run_again_anywhere , x=lossmargin, c=0, level=90))


sler_next4_female_cct <- with(subset(sler4,female==1), 
                        rdrobust(y=run_again_next4yrs , x=lossmargin, c=0))
summary(sler_next4_female_cct) # female cands, effect of 0.336
sler_next4_female_cct_90 <- with(subset(sler4,female==1), 
                           rdrobust(y=run_again_next4yrs , x=lossmargin, c=0, level=90))

sler_next4_male_cct <- with(subset(sler4,female==0), 
                      rdrobust(y=run_again_next4yrs , x=lossmargin, c=0))
summary(sler_next4_male_cct) # male cands, effect of 0.304
sler_next4_male_cct_90 <- with(subset(sler4,female==0), 
                         rdrobust(y=run_again_next4yrs , x=lossmargin, c=0, level=90))



# significance test of difference in coefs:
(sler_diff <- sler_female_cct$coef[1,1] - sler_male_cct$coef[1,1]) # 0.031
Z_gender = (sler_female_cct$coef[1,1] - sler_male_cct$coef[1,1])/
	sqrt(sler_female_cct$se[3,1]^2 + sler_male_cct$se[3,1]^2)
(Pvalue_gender <- 2*pnorm(-abs(Z_gender))) # 0.18


sler_firsttime_female_cct <- with(subset(sler_firsttime,female==1), 
                                        rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
sler_firsttime_female_cct_90 <- with(subset(sler_firsttime,female==1), 
                                           rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(sler_firsttime_female_cct) # female cands, effect of 0.477

sler_firsttime_male_cct <- with(subset(sler_firsttime,female==0), 
                                      rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
sler_firsttime_male_cct_90 <- with(subset(sler_firsttime,female==0), 
                                         rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(sler_firsttime_male_cct) # male cands, effect of 0.453



sler_next4_firsttime_female_cct <- with(subset(sler_firsttime,female==1), 
                                        rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
sler_next4_firsttime_female_cct_90 <- with(subset(sler_firsttime,female==1), 
                                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(sler_next4_firsttime_female_cct) # female cands, effect of 0.493

sler_next4_firsttime_male_cct <- with(subset(sler_firsttime,female==0), 
                                      rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
sler_next4_firsttime_male_cct_90 <- with(subset(sler_firsttime,female==0), 
                                         rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(sler_next4_firsttime_male_cct) # male cands, effect of 0.480



# significance test of difference in coefs:
(sler_diff <- sler_next4_firsttime_female_cct$coef[1,1] - sler_next4_firsttime_male_cct$coef[1,1]) # 0.013
Z_gender = (sler_next4_firsttime_female_cct$coef[1,1] - sler_next4_firsttime_male_cct$coef[1,1])/
  sqrt(sler_next4_firsttime_female_cct$se[3,1]^2 + sler_next4_firsttime_male_cct$se[3,1]^2)
(Pvalue_gender_sler <- 2*pnorm(-abs(Z_gender))) # 0.633



#---------------------#
#### CEDA analysis ####
#---------------------#
ceda_all_cct <- with(subset(ceda4), 
										 rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname))
ceda_all_cct_90 <- with(subset(ceda4), 
                     rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_all_cct) # across all cands, effect of 0.167


## next 4 years:
ceda_next4_all_cct <- with(subset(ceda4), 
                     rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_all_cct_90 <- with(subset(ceda4), 
                        rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))

## next 4 years, first-time only:
ceda_next4_firsttime_all_cct <- with(subset(ceda_firsttime), 
                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_firsttime_all_cct_90 <- with(subset(ceda_firsttime), 
                              rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))


## by gender:
ceda_female_cct <- with(subset(ceda4,female==1), 
												rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname))
ceda_female_cct_90 <- with(subset(ceda4,female==1), 
                        rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_female_cct) # female cands, effect of 0.205


ceda_male_cct <- with(subset(ceda4,female==0), 
											rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname))
ceda_male_cct_90 <- with(subset(ceda4,female==0), 
                      rdrobust(y=run_again_anywhere, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_male_cct) # male cands, effect of 0.169

# significance test of difference in coefs:
(ceda_diff <- ceda_female_cct$coef[1,1] - ceda_male_cct$coef[1,1]) # 0.013
Z_gender = (ceda_female_cct$coef[1,1] - ceda_male_cct$coef[1,1])/
  sqrt(ceda_female_cct$se[3,1]^2 + ceda_male_cct$se[3,1]^2)
(Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))) # 0.745


ceda_next4_female_cct <- with(subset(ceda4,female==1), 
                        rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_female_cct_90 <- with(subset(ceda4,female==1), 
                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_next4_female_cct) # female cands, effect of 0.1

ceda_next4_male_cct <- with(subset(ceda4,female==0), 
                      rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_male_cct_90 <- with(subset(ceda4,female==0), 
                         rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_next4_male_cct) # male cands, effect of 0.11

(ceda_diff <- ceda_next4_female_cct$coef[1,1] - ceda_next4_male_cct$coef[1,1]) # 0.008
Z_gender = (ceda_next4_female_cct$coef[1,1] - ceda_next4_male_cct$coef[1,1])/
  sqrt(ceda_next4_female_cct$se[3,1]^2 + ceda_next4_male_cct$se[3,1]^2)
(Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))) # 0.825



#### CEDA analysis, robustness checks ####
## Using FEs and alternative bandwidths other than those from rdrobust:
ceda_runanywhere_female_lm_wasserman <- felm(run_again_anywhere ~ lossmargin*lost | 0 | 0 | cntyname,
                                   data=filter(ceda4, female==1 & lossmargin>=(-0.10) & lossmargin<=(+0.10)))
summary(ceda_runanywhere_female_lm_wasserman) # 0.171

ceda_runanywhere_male_lm_wasserman <- felm(run_again_anywhere ~ lossmargin*lost | 0 | 0 | cntyname,
                                 data=filter(ceda4, female==0 & lossmargin>=(-0.094) & lossmargin<=(+0.094)))
summary(ceda_runanywhere_male_lm_wasserman) # 0.128



ceda_runanywhere_wasserman_nofes <- felm(run_again_anywhere ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                                   data=filter(ceda4, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

ceda_runanywhere_wasserman <- felm(run_again_anywhere ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                             data=filter(ceda4, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

summary(ceda_runanywhere_wasserman) # bigger effect for female cands: intx of lost*female of -0.04

ceda_next4_wasserman <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                             data=filter(ceda4, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

summary(ceda_next4_wasserman) # bigger effect for female cands: insig intx of -0.03

ceda_next4_firsttime_wasserman <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                             data=filter(ceda_firsttime, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

summary(ceda_next4_firsttime_wasserman) # bigger effect for female cands: insig intx of -0.04

ceda_next4_wasserman_nofe <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                             data=filter(ceda4, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

summary(ceda_next4_wasserman_nofe) # bigger effect for female cands: insig intx of -0.03

ceda_next4_firsttime_wasserman_nofe <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                                  data=filter(ceda_firsttime, lossmargin>=(-0.096) & lossmargin<=(+0.096)))

summary(ceda_next4_firsttime_wasserman_nofe) # bigger effect for female cands: insig intx of -0.03


ceda_runanywhere_female_lm_optimal <- felm(run_again_anywhere ~ lossmargin*lost | 0 | 0 | cntyname,
                                           data=filter(ceda4, female==1 & lossmargin>=(-ceda_female_cct$bws[1,1]) & lossmargin<=(+ceda_female_cct$bws[1,2])))
summary(ceda_runanywhere_female_lm_optimal) # 0.125

ceda_runanywhere_male_lm_optimal <- felm(run_again_anywhere ~ lossmargin*lost | 0 | 0 | cntyname,
                                         data=filter(ceda4, female==0 & lossmargin>=(-ceda_male_cct$bws[1,1]) & lossmargin<=(+ceda_male_cct$bws[1,2])))
summary(ceda_runanywhere_male_lm_optimal) # 0.134



ceda_runanywhere_optimal_nofes <- felm(run_again_anywhere ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                                       data=filter(ceda4, lossmargin>=(-ceda_all_cct$bws[1,1]) & lossmargin<=(+ceda_all_cct$bws[1,2])))

ceda_runanywhere_optimal <- felm(run_again_anywhere ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                                 data=filter(ceda4, lossmargin>=(-ceda_all_cct$bws[1,1]) & lossmargin<=(+ceda_all_cct$bws[1,2])))

summary(ceda_runanywhere_optimal) # insig intx of lost*female of -0.007


ceda_next4_optimal_nofe <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                                data=filter(ceda4, lossmargin>=(-ceda_next4_all_cct$bws[1,1]) & lossmargin<=(+ceda_next4_all_cct$bws[1,2])))
summary(ceda_next4_optimal_nofe) # insig intx of 0.005

ceda_next4_optimal <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                           data=filter(ceda4, lossmargin>=(-ceda_next4_all_cct$bws[1,1]) & lossmargin<=(+ceda_next4_all_cct$bws[1,2])))
summary(ceda_next4_optimal) # insig intx of -0.001

ceda_next4_firsttime_optimal <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | cntyname + year | 0 | cntyname, 
                                     data=filter(ceda_firsttime, lossmargin>=(-ceda_next4_firsttime_all_cct$bws[1,1]) & lossmargin<=(+ceda_next4_firsttime_all_cct$bws[1,2])))
summary(ceda_next4_firsttime_optimal) # bigger effect for female cands: insig intx of -0.039

ceda_next4_firsttime_optimal_nofe <- felm(run_again_next4yrs ~ lossmargin*(lost)*female | 0 | 0 | cntyname, 
                                          data=filter(ceda_firsttime, lossmargin>=(-ceda_next4_firsttime_all_cct$bws[1,1]) & lossmargin<=(+ceda_next4_firsttime_all_cct$bws[1,2])))
summary(ceda_next4_firsttime_optimal_nofe) # insig intx of -0.035



### rdrobust versions:
ceda_next4_firsttime_female_cct <- with(subset(ceda_firsttime,female==1), 
                              rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_firsttime_female_cct_90 <- with(subset(ceda_firsttime,female==1), 
                                 rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_next4_firsttime_female_cct) # female cands, effect of 0.096

ceda_next4_firsttime_male_cct <- with(subset(ceda_firsttime,female==0), 
                            rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname))
ceda_next4_firsttime_male_cct_90 <- with(subset(ceda_firsttime,female==0), 
                               rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, cluster=cntyname, level=90))
summary(ceda_next4_firsttime_male_cct) # male cands, effect of 0.074


(ceda_diff <- ceda_next4_firsttime_female_cct$coef[1,1] - ceda_next4_firsttime_male_cct$coef[1,1]) # 0.021
Z_gender = (ceda_next4_firsttime_female_cct$coef[1,1] - ceda_next4_firsttime_male_cct$coef[1,1])/
  sqrt(ceda_next4_firsttime_female_cct$se[3,1]^2 + ceda_next4_firsttime_male_cct$se[3,1]^2)
(Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))) # 0.572


stargazer(ceda_runanywhere_female_lm_wasserman,ceda_runanywhere_male_lm_wasserman,
          ceda_runanywhere_wasserman_nofes,ceda_runanywhere_wasserman,
          ceda_next4_wasserman_nofe,ceda_next4_wasserman,
          ceda_next4_firsttime_wasserman_nofe,ceda_next4_firsttime_wasserman,
          out = "tables/wasserman.tex",
          dep.var.labels = c("Run again anytime","Run again in next 4 years"),
          column.labels = c("Women Only","Men Only","All candidates","First-time candidates"),
          column.separate = c(1,1,4,2),
          omit.stat = c("ser","f","adj.rsq"),
          omit = "Constant",
          covariate.labels = c("Loss margin","Lost","Female","Loss margin $\\times$ lost","Loss margin $\\times$ Female","Lost $\\times$ Female","Loss margin $\\times$ Lost $\\times$ Female"),
          add.lines = list(c("County+year fixed effects?","No","No","No","Yes","No","Yes","No","Yes"),
                           c("Bandwidth (Wasserman's)",0.10, 0.094, 0.096, 0.096, 0.096, 0.096, 0.096, 0.096, 0.096)),
          notes="Standard errors clustered by county and uniform kernel used for local linear estimation in all models.",
          float = F)


stargazer(ceda_runanywhere_female_lm_optimal,ceda_runanywhere_male_lm_optimal,
          ceda_runanywhere_optimal_nofes,ceda_runanywhere_optimal,
          ceda_next4_optimal_nofe,ceda_next4_optimal,
          ceda_next4_firsttime_optimal_nofe,ceda_next4_firsttime_optimal,
          out = "tables/wasserman_optimalbw.tex",
          dep.var.labels = c("Run again anytime","Run again in next 4 years"),
          column.labels = c("Women Only","Men Only","All candidates","First-time candidates"),
          column.separate = c(1,1,4,2),
          omit.stat = c("ser","f","adj.rsq"),
          omit = "Constant",
          covariate.labels = c("Loss margin","Lost","Female","Loss margin $\\times$ lost","Loss margin $\\times$ Female","Lost $\\times$ Female","Loss margin $\\times$ Lost $\\times$ Female"),
          add.lines = list(c("County+year fixed effects?","No","No","No","Yes","No","Yes","No","Yes"),
                           c("Bandwidth (MSE-Optimal)",
                             round(ceda_female_cct$bws[1, 1],3), round(ceda_male_cct$bws[1, 1],3),
                             round(ceda_all_cct$bws[1, 1],3), round(ceda_all_cct$bws[1, 1],3),
                             round(ceda_next4_all_cct$bws[1, 1],3), round(ceda_next4_all_cct$bws[1, 1],3),
                             round(ceda_next4_firsttime_all_cct$bws[1, 1],3), round(ceda_next4_firsttime_all_cct$bws[1, 1],3))),
          notes="Standard errors clustered by county and uniform kernel used for local linear estimation in all models.",
          float = F)


rdd_ceda_tab <- rd.export4(list(
  ceda_male_cct,
  ceda_female_cct,
  
  ceda_firsttime_male_cct,
  ceda_firsttime_female_cct,
  
  ceda_next4_male_cct,
  ceda_next4_female_cct,
  
  ceda_next4_firsttime_male_cct,
  ceda_next4_firsttime_female_cct
)) %>% as.data.frame()

rdd_ceda_tab$Subset <- c("Male","","Female","","Male, first-time","","Female, first-time","","Male","","Female","","Male, first-time","","Female, first-time","")
rdd_ceda_tab$DV <- c("Run Anytime","","Run Anytime","","Run Anytime","","Run Anytime","","Run in next 4 years","","Run in next 4 years","","Run in next 4 years","","Run in next 4 years","")
rdd_ceda_tab <- select(rdd_ceda_tab,
                  DV,Subset,Coef,`p-value`,BW,Obs)

print(xtable(rdd_ceda_tab),file = "tables/rdd_ceda_bygender_wasserman.tex",floating = F,include.rownames = F)


#### CEDA: alternative bandwidths ####

coefs.female <- NULL
coefs.male <- NULL
diffs_runanytime <- NULL
coefs.next4.female <- NULL
coefs.next4.male <- NULL
diffs_next4 <- NULL
coefs.next4.firsttime.female <- NULL
coefs.next4.firsttime.male <- NULL
diffs_next4_firsttime <- NULL
for(i in seq(0.01,0.5,0.01)){
  
  fit_female_cct_temp <- with(filter(ceda4,female==1),
                                                   rdrobust(y = run_again_anywhere,x = lossmargin,
                                                            c = 0,cluster = cntyname,
                                                            h = i
                                                   ))
  coefs.female <- bind_rows(coefs.female,rd.export.numeric(list(fit_female_cct_temp)))
  
  fit_male_cct_temp <- with(filter(ceda4,female==0),
                                                      rdrobust(y = run_again_anywhere,x = lossmargin,
                                                               c = 0,
                                                               cluster = cntyname,
                                                               h = i
                                                      ))
  coefs.male <- bind_rows(coefs.male,rd.export.numeric(list(fit_male_cct_temp)))
  
  ceda_diff <- fit_female_cct_temp$coef[1,1] - fit_male_cct_temp$coef[1,1]
  Z_gender = (fit_female_cct_temp$coef[1,1] - fit_male_cct_temp$coef[1,1])/
    sqrt(fit_female_cct_temp$se[3,1]^2 + fit_male_cct_temp$se[3,1]^2)
  Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))
  diffs_runanytime <- bind_rows(diffs_runanytime,bind_cols(ceda_diff,Pvalue_gender_ceda))
  
  
  fit_next4_female_cct_temp <- with(filter(ceda4,female==1),
                              rdrobust(y = run_again_next4yrs,x = lossmargin,
                                       c = 0,cluster = cntyname,
                                       h = i
                              ))
  coefs.next4.female <- bind_rows(coefs.next4.female,rd.export.numeric(list(fit_next4_female_cct_temp)))
  
  fit_next4_male_cct_temp <- with(filter(ceda4,female==0),
                            rdrobust(y = run_again_next4yrs,x = lossmargin,
                                     c = 0,
                                     cluster = cntyname,
                                     h = i
                            ))
  coefs.next4.male <- bind_rows(coefs.next4.male,rd.export.numeric(list(fit_next4_male_cct_temp)))
  
  ceda_diff <- fit_next4_female_cct_temp$coef[1,1] - fit_next4_male_cct_temp$coef[1,1]
  Z_gender = (fit_next4_female_cct_temp$coef[1,1] - fit_next4_male_cct_temp$coef[1,1])/
    sqrt(fit_next4_female_cct_temp$se[3,1]^2 + fit_next4_male_cct_temp$se[3,1]^2)
  Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))
  diffs_next4 <- bind_rows(diffs_next4,bind_cols(ceda_diff,Pvalue_gender_ceda))
  
  fit_next4_firsttime_female_cct_temp <- with(filter(ceda_firsttime,female==1),
                                    rdrobust(y = run_again_next4yrs,x = lossmargin,
                                             c = 0,cluster = cntyname,
                                             h = i
                                    ))
  coefs.next4.firsttime.female <- bind_rows(coefs.next4.firsttime.female,rd.export.numeric(list(fit_next4_firsttime_female_cct_temp)))
  
  fit_next4_firsttime_male_cct_temp <- with(filter(ceda_firsttime,female==0),
                                  rdrobust(y = run_again_next4yrs,x = lossmargin,
                                           c = 0,
                                           cluster = cntyname,
                                           h = i
                                  ))
  coefs.next4.firsttime.male <- bind_rows(coefs.next4.firsttime.male,rd.export.numeric(list(fit_next4_firsttime_male_cct_temp)))
  
  ceda_diff <- fit_next4_firsttime_female_cct_temp$coef[1,1] - fit_next4_firsttime_male_cct_temp$coef[1,1]
  Z_gender = (fit_next4_firsttime_female_cct_temp$coef[1,1] - fit_next4_firsttime_male_cct_temp$coef[1,1])/
    sqrt(fit_next4_firsttime_female_cct_temp$se[3,1]^2 + fit_next4_firsttime_male_cct_temp$se[3,1]^2)
  Pvalue_gender_ceda <- 2*pnorm(-abs(Z_gender))
  diffs_next4_firsttime <- bind_rows(diffs_next4_firsttime,bind_cols(ceda_diff,Pvalue_gender_ceda))
  
}

coefs.female$gender <- "Female"
coefs.male$gender <- "Male"
coefs.both <- bind_rows(coefs.female,coefs.male)
coefs.next4.female$gender <- "Female"
coefs.next4.male$gender <- "Male"
coefs.next4.both <- bind_rows(coefs.next4.female,coefs.next4.male)
coefs.next4.firsttime.female$gender <- "Female"
coefs.next4.firsttime.male$gender <- "Male"
coefs.next4.firsttime.both <- bind_rows(coefs.next4.firsttime.female,coefs.next4.firsttime.male)


colnames(diffs_runanytime)[1:2] <- c("Gender diff.","p-value of difference")
colnames(diffs_next4)[1:2] <- c("Gender diff.","p-value of difference")
colnames(diffs_next4_firsttime)[1:2] <- c("Gender diff.","p-value of difference")
diffs_runanytime$bw <- seq(0.01,0.5,0.01)
diffs_next4$bw <- seq(0.01,0.5,0.01)
diffs_next4_firsttime$bw <- seq(0.01,0.5,0.01)

pdf("figures/bandwidths_runanytime.pdf",width=7,height=5)
ggplot(data = coefs.both) +
  geom_point(aes(x=bw,y=coef,color=gender),size = 3)+
  geom_errorbar(aes(x = bw,color=gender, 
                    ymin = cilo, ymax = cihi, height=.1)) +
  geom_text(data=filter(diffs_runanytime,`p-value of difference`<0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="red",size=1.5) +
  geom_text(data=filter(diffs_runanytime,`p-value of difference`>=0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="black",size=1.5) + 
  scale_x_continuous("Bandwidth",breaks=seq(0,50,10)) + 
  ylab("Effect of loss on running\nagain anytime") +
  scale_color_manual("Candidate gender",breaks = c("Female","Male"),values=c("#440154FF","#21908CFF")) + 
  coord_cartesian(ylim=c(-0.35,0.08),xlim=c(0,55)) + 
  geom_hline(yintercept = 0,size=.5,colour="red",linetype="dotted") +
  theme_bw() +
  theme(axis.text = element_text(size=15),axis.title = element_text(size=15),
        axis.text.y = element_text(angle = 0, hjust = 0, color="black"),
        legend.position = c(0.6,0.1))
dev.off()



pdf("figures/bandwidths_next4.pdf",width=7,height=5)
ggplot(data = coefs.next4.both) +
  geom_point(aes(x=bw,y=coef,color=gender),size = 3)+
  geom_errorbar(aes(x = bw,color=gender, 
                    ymin = cilo, ymax = cihi, height=.1)) +
  geom_text(data=filter(diffs_next4,`p-value of difference`<0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="red",size=1.5) +
  geom_text(data=filter(diffs_next4,`p-value of difference`>=0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="black",size=1.5) + 
  scale_x_continuous("Bandwidth",breaks=seq(0,50,10)) + 
  ylab("Effect of loss on running\nagain in next 4 years") +
  scale_color_manual("Candidate gender",breaks = c("Female","Male"),values=c("#440154FF","#21908CFF")) + 
  coord_cartesian(ylim=c(-0.35,0.07),xlim=c(0,55)) + 
  geom_hline(yintercept = 0,size=.5,colour="red",linetype="dotted") +
  theme_bw() +
  theme(axis.text = element_text(size=15),axis.title = element_text(size=15),
        axis.text.y = element_text(angle = 0, hjust = 0, color="black"),
        legend.position = c(0.6,0.1))
dev.off()


pdf("figures/bandwidths_next4_firsttime.pdf",width=7,height=5)
ggplot(data = coefs.next4.firsttime.both) +
  geom_point(aes(x=bw,y=coef,color=gender),size = 3)+
  geom_errorbar(aes(x = bw,color=gender, 
                    ymin = cilo, ymax = cihi, height=.1)) +
  # geom_text(data=filter(diffs_next4_firsttime,`p-value of difference`<0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="red",size=1.5) + 
  geom_text(data=filter(diffs_next4_firsttime,`p-value of difference`>=0.1),aes(x=bw*100,label=paste0("Diff = ",round(`Gender diff.`,2),", p-value of diff = ",round(`p-value of difference`,3))),y=0.01,angle=45,hjust=0,col="black",size=1.5) + 
  scale_x_continuous("Bandwidth",breaks=seq(0,50,10)) + 
  ylab("Effect of loss on running again in\nnext 4 years (first-time candidates)") +
  scale_color_manual("Candidate gender",breaks = c("Female","Male"),values=c("#440154FF","#21908CFF")) + 
  coord_cartesian(ylim=c(-0.35,0.07),xlim=c(0,55)) + 
  geom_hline(yintercept = 0,size=.5,colour="red",linetype="dotted") +
  theme_bw() +
  theme(axis.text = element_text(size=15),axis.title = element_text(size=15),
        axis.text.y = element_text(angle = 0, hjust = 0, color="black"),
        legend.position = c(0.6,0.1))
dev.off()





#------------------------#
#### Mayoral analysis ####
#------------------------#

mayors_all_cct <- with(subset(mayors4), 
											 rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
mayors_all_cct_90 <- with(subset(mayors4), 
                       rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(mayors_all_cct) # effect of 0.091


mayors_next4_all_cct <- with(subset(mayors4), 
                       rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
mayors_next4_all_cct_90 <- with(subset(mayors4), 
                          rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))


mayors_female_cct <- with(subset(mayors4, female==1), 
													rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
mayors_female_cct_90 <- with(subset(mayors4, female==1), 
                          rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(mayors_female_cct) # effect of 0.319

mayors_male_cct <- with(subset(mayors4, female==0), 
												rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
mayors_male_cct_90 <- with(subset(mayors4, female==0), 
                        rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(mayors_male_cct) # effect of 0.198


mayors_next4_female_cct <- with(subset(mayors4, female==1), 
                          rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
mayors_next4_female_cct_90 <- with(subset(mayors4, female==1), 
                             rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(mayors_next4_female_cct) # effect of 0.305

mayors_next4_male_cct <- with(subset(mayors4, female==0), 
                        rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
mayors_next4_male_cct_90 <- with(subset(mayors4, female==0), 
                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(mayors_next4_male_cct) # effect of 0.219


# significance test of difference in coefs:
(mayors_diff <- mayors_female_cct$coef[1,1] - mayors_male_cct$coef[1,1]) # 0.070
Z_gender = (mayors_female_cct$coef[1,1] - mayors_male_cct$coef[1,1])/
	sqrt(mayors_female_cct$se[3,1]^2 + mayors_male_cct$se[3,1]^2)
(Pvalue_gender_mayors <- 2*pnorm(-abs(Z_gender))) # 0.117


(mayors_next4_female_cct$coef[1,1] - mayors_next4_male_cct$coef[1,1]) # 0.086
Z_gender = (mayors_next4_female_cct$coef[1,1] - mayors_next4_male_cct$coef[1,1])/
  sqrt(mayors_next4_female_cct$se[3,1]^2 + mayors_next4_male_cct$se[3,1]^2)
(Pvalue_gender <- 2*pnorm(-abs(Z_gender))) # 0.281


mayors_firsttime_female_cct <- with(subset(mayors_firsttime,female==1), 
                                          rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
mayors_firsttime_female_cct_90 <- with(subset(mayors_firsttime,female==1), 
                                             rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(mayors_firsttime_female_cct) # female cands, effect of 0.429

mayors_firsttime_male_cct <- with(subset(mayors_firsttime,female==0), 
                                        rdrobust(y=run_again_anywhere, x=lossmargin, c=0))
mayors_firsttime_male_cct_90 <- with(subset(mayors_firsttime,female==0), 
                                           rdrobust(y=run_again_anywhere, x=lossmargin, c=0, level=90))
summary(mayors_firsttime_male_cct) # male cands, effect of 0.355


mayors_next4_firsttime_female_cct <- with(subset(mayors_firsttime,female==1), 
                                          rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
mayors_next4_firsttime_female_cct_90 <- with(subset(mayors_firsttime,female==1), 
                                             rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(mayors_next4_firsttime_female_cct) # female cands, effect of 0.435

mayors_next4_firsttime_male_cct <- with(subset(mayors_firsttime,female==0), 
                                        rdrobust(y=run_again_next4yrs, x=lossmargin, c=0))
mayors_next4_firsttime_male_cct_90 <- with(subset(mayors_firsttime,female==0), 
                                           rdrobust(y=run_again_next4yrs, x=lossmargin, c=0, level=90))
summary(mayors_next4_firsttime_male_cct) # male cands, effect of 0.381



# significance test of difference in coefs:
(mayors_diff <- mayors_next4_firsttime_female_cct$coef[1,1] - mayors_next4_firsttime_male_cct$coef[1,1]) # 0.054
Z_gender = (mayors_next4_firsttime_female_cct$coef[1,1] - mayors_next4_firsttime_male_cct$coef[1,1])/
  sqrt(mayors_next4_firsttime_female_cct$se[3,1]^2 + mayors_next4_firsttime_male_cct$se[3,1]^2)
(Pvalue_gender_mayors <- 2*pnorm(-abs(Z_gender))) # 0.505



# ------------------- #
#### Figures: SLER ####
# ------------------- #

### SLERs - descriptive
sler_summ <- sler4 %>%
  group_by(voteshare>=0.5) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(lose = !win)

plot_descriptive_sler <- ggplot(sler_summ, aes(x=lose,y=run_again_anywhere)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.04,label=round(run_again_anywhere*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal()
ggsave(plot_descriptive_sler,filename = "figures/descriptive_sler.pdf",width=3,height=3)


plot_descriptive3_sler <- ggplot(sler_summ, aes(x=lose,y=run_again_next4yrs)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.04,label=round(run_again_next4yrs*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal()
ggsave(plot_descriptive3_sler,filename = "figures/descriptive3_sler.pdf",width=3,height=3)

# prep for binned plot:
width <- .01
data.graph <- sler4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere)))

plot_descriptive_cont_sler <- ggplot(data.graph.all) + 
  geom_point(aes(x=mid,y=bin.mean,size=n),position=position_jitter(height=0.05),shape=1) + 
  geom_smooth(data=filter(sler4,voteshare<0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  geom_smooth(data=filter(sler4,voteshare>=0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal() + 
  theme(legend.position = "bottom")
ggsave(plot_descriptive_cont_sler,filename = "figures/descriptive_cont_sler.pdf",width=4,height=4)

# with gender:
width <- .01
data.graph <- sler4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere,female) %>%
  group_by(bin,mid,female) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere)))


sler4 <- sler4 %>%
  mutate(female2=factor(female,levels = c(FALSE,TRUE),labels=c("Male","Female")))

plot_descriptive_cont_nopts_nose_gender_sler <- ggplot(data.graph.all) + 
  # geom_point(data=filter(data.graph.all,female==1),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#440154FF",alpha=0.8) + 
  # geom_point(data=filter(data.graph.all,female==0),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#21908CFF",alpha=0.8) + 
  geom_smooth(data=filter(sler4,voteshare<0.5 & female2=="Female"),aes(x=voteshare,y=run_again_anywhere,color=female2),method="loess",se=F) + 
  geom_smooth(data=filter(sler4,voteshare>=0.5 & female2=="Female"),aes(x=voteshare,y=run_again_anywhere,color=female2),method="loess",se=F) + 
  geom_smooth(data=filter(sler4,voteshare<0.5 & female2=="Male"),aes(x=voteshare,y=run_again_anywhere,color=female2),method="loess",lty=2,se=F) + 
  geom_smooth(data=filter(sler4,voteshare>=0.5 & female2=="Male"),aes(x=voteshare,y=run_again_anywhere,color=female2),method="loess",lty=2,se=F) + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_color_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal() + 
  theme(legend.position = c(0.8,0.15))
ggsave(plot_descriptive_cont_nopts_nose_gender_sler,filename = "figures/descriptive_cont_nopts_nose_gender_sler.pdf",width=4,height=4)



sler_summ_gender <- sler4 %>%
  group_by(voteshare>=0.5,female) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(female=factor(female,levels = c(FALSE,TRUE),labels=c("Male","Female")),
         lose = !win) %>%
  ungroup()


plot_descriptive_sler_gender <- ggplot(sler_summ_gender, aes(x=lose,y=run_again_anywhere,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.03,label=round(run_again_anywhere*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive_sler_gender,filename = "figures/descriptive_gender_sler.pdf",width=6,height=4)


plot_descriptive3_sler_gender <- ggplot(sler_summ_gender, aes(x=lose,y=run_again_next4yrs,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.03,label=round(run_again_next4yrs*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("Nationwide State Legislative") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive3_sler_gender,filename = "figures/descriptive3_gender_sler.pdf",width=6,height=4)



# prep data for binned means:
width <- .005
data.graph <- sler4 %>%
	mutate(bin=cut(lossmargin, breaks=seq(-0.5,0.5, width))
	)
bins <- data.frame(bin = levels(data.graph$bin),
									 mid = seq(-0.5 + width/2, 0.5 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.female <- data.graph %>% 
	filter(female==1) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.male <- data.graph %>% 
	filter(female==0) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin, mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))


pdf("figures/rdplot_runanywhere_all_slers_binned.pdf", width=3, height=3.5)
(ggplot(sler4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# all candidates:
		geom_smooth(data=subset(sler4, voteshare_toptwo > 0.5),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_smooth(data=subset(sler4, voteshare_toptwo < 0.5),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_point(data=data.graph.all, 
							 aes(x=mid, y=bin.mean, size=n), pch=1, color="black") + 
		# geom_vline(xintercept=0.5) + 
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013)) +
		scale_x_continuous(breaks=seq(-0.40, 0.40, 0.02),
											 limits=c(-sler_all_cct$bws["h", "left"],
											 				 sler_all_cct$bws["h", "right"]),
											 labels=seq(0.4, -0.4, -0.020)) + 
		theme_minimal() + 
		labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("Nationwide State Legislative") +
    theme(plot.title = element_text(hjust = 0),
          legend.position = "bottom", legend.text = element_text(size=4.5),legend.title = element_text(size=4.5),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    annotate('text', x = (0.465-0.5), y = .65, 0.05, hjust=0, parse=TRUE, color="black",
						 label = paste('hat(tau)==', round(sler_all_cct$coef['Conventional', ], 2),
						 							'~(list(', round(sler_all_cct$ci['Robust', 1], 2),
						 							',', round(sler_all_cct$ci['Robust', 2], 2), '))'))
)
dev.off()


pdf("figures/rdplot_runanywhere_gender_slers_binned.pdf", width=3, height=3.5)
(ggplot(sler4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# female candidates:
		geom_smooth(data=subset(sler4, voteshare_toptwo > 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_smooth(data=subset(sler4, voteshare_toptwo < 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_point(data=data.graph.female, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#440154FF") + 
		# male candidates:
		geom_smooth(data=subset(sler4, voteshare_toptwo > 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_smooth(data=subset(sler4, voteshare_toptwo < 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_point(data=data.graph.male, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#21908CFF") + 
		# geom_vline(xintercept=0.5) + 
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013),xlim=c(-sler_female_cct$bws["h", "left"],
		                                              sler_female_cct$bws["h", "right"])) +
		scale_x_continuous(breaks=seq(-0.4, 0.4, .02),
											 limits=c(-sler_female_cct$bws["h", "left"]-0.01,
											 				 sler_female_cct$bws["h", "right"]+0.01),
											 labels=seq(0.4, -0.4, -.02)) + 
    theme_minimal() + 
    theme(plot.title = element_text(hjust = 0),
		      legend.position = "bottom", legend.text = element_text(size=5),legend.title = element_text(size=7.5),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("Nationwide State Legislative") +
		annotate('text', x = -0.054, y = .6, 0.05, hjust=0, parse=TRUE, color="#21908CFF",
						 label = paste('Male:~hat(tau)==', round(sler_male_cct$coef['Conventional', ], 2),
						 							'~(list(', round(sler_male_cct$ci['Robust', 1], 2),
						 							',', round(sler_male_cct$ci['Robust', 2], 2), '))')) + 
		annotate('text', x = -0.054, y = .7, 0.05, hjust=0, parse=TRUE, color="#440154FF",
						 label = paste('Female:~hat(tau)==', round(sler_female_cct$coef['Conventional', ], 2),
						 							'~(list(', round(sler_female_cct$ci['Robust', 1], 2),
						 							',', round(sler_female_cct$ci['Robust', 2], 2), '))')) 
)
dev.off()



# ------------------- #
#### Figures: CEDA ####
# ------------------- #

### CEDA - descriptive
ceda_summ <- ceda4 %>%
  group_by(voteshare>=0.5) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(lose = !win)

plot_descriptive_ceda <- ggplot(ceda_summ, aes(x=lose,y=run_again_anywhere)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.04,label=round(run_again_anywhere*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("CA Local") +
  theme_minimal()
ggsave(plot_descriptive_ceda,filename = "figures/descriptive_ceda.pdf",width=3,height=3)


plot_descriptive3_ceda <- ggplot(ceda_summ, aes(x=lose,y=run_again_next4yrs)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.04,label=round(run_again_next4yrs*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("CA Local") +
  theme_minimal()
ggsave(plot_descriptive3_ceda,filename = "figures/descriptive3_ceda.pdf",width=3,height=3)

# prep for binned plot:
width <- .01
data.graph <- ceda4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere)))

plot_descriptive_cont_ceda <- ggplot(data.graph.all) + 
  geom_point(aes(x=mid,y=bin.mean,size=n),position=position_jitter(height=0.05),shape=1) + 
  geom_smooth(data=filter(ceda4,voteshare<0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  geom_smooth(data=filter(ceda4,voteshare>=0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("CA Local") +
  theme_minimal() + 
  theme(legend.position = "bottom")
ggsave(plot_descriptive_cont_ceda,filename = "figures/descriptive_cont_ceda.pdf",width=4,height=4)


## continuous descriptive, by gender:
width <- .01
data.graph <- ceda4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere,female) %>%
  group_by(bin,mid,female) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere))) %>%
  na.omit() %>%
  mutate(female=factor(female,levels = c(0,1),labels=c("Male","Female")))


ceda4 <- ceda4 %>%
  mutate(female2=factor(female,levels = c(0,1),labels=c("Male","Female")))

plot_descriptive_cont_nopts_nose_gender_ceda <- ggplot(data.graph.all) + 
  # geom_point(data=filter(data.graph.all,female=="Male"),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#21908CFF",alpha=0.8) + 
  # geom_point(data=filter(data.graph.all,female=="Female"),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#440154FF",alpha=0.8) + 
  geom_smooth(data=filter(ceda4,voteshare<0.5 & female==1),aes(x=voteshare,y=run_again_anywhere,color=female2),
              method="loess",se=F) + 
  geom_smooth(data=filter(ceda4,voteshare>=0.5 & female==1),aes(x=voteshare,y=run_again_anywhere,color=female2),
              method="loess",se=F) + 
  geom_smooth(data=filter(ceda4,voteshare<0.5 & female==0),aes(x=voteshare,y=run_again_anywhere,color=female2),
              method="loess",lty=2,se=F) + 
  geom_smooth(data=filter(ceda4,voteshare>=0.5 & female==0),aes(x=voteshare,y=run_again_anywhere,color=female2),
              method="loess",lty=2,se=F) + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_color_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("CA Local") +
  theme_minimal() + 
  theme(legend.position = c(0.8,0.15))
ggsave(plot_descriptive_cont_nopts_nose_gender_ceda,filename = "figures/descriptive_cont_nopts_nose_gender_ceda.pdf",width=4,height=4)


ceda_summ_gender <- ceda4 %>%
  group_by(voteshare>=0.5,female) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(female=factor(female,levels = c(0,1),labels=c("Male","Female")),
         lose=!win) %>%
  ungroup()


plot_descriptive_ceda_gender <- ggplot(ceda_summ_gender, aes(x=lose,y=run_again_anywhere,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.03,label=round(run_again_anywhere*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("CA Local") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive_ceda_gender,filename = "figures/descriptive_gender_ceda.pdf",width=6,height=4)


plot_descriptive3_ceda_gender <- ggplot(ceda_summ_gender, aes(x=lose,y=run_again_next4yrs,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.03,label=round(run_again_next4yrs*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("CA Local") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive3_ceda_gender,filename = "figures/descriptive3_gender_ceda.pdf",width=6,height=4)






# prep data for binned means:
width <- .005
data.graph <- ceda4 %>%
	mutate(bin=cut(lossmargin, breaks=seq(-0.5,0.5, width))
	)
bins <- data.frame(bin = levels(data.graph$bin),
									 mid = seq(-0.5 + width/2, 0.5 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.female <- data.graph %>% 
	filter(female==1) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.male <- data.graph %>% 
	filter(female==0) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin, mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

pdf("figures/rdplot_runanywhere_all_ceda_binned.pdf", width=3, height=3.5)
(ggplot(ceda4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# all candidates:
		geom_smooth(data=subset(ceda4, voteshare_toptwo > 0.5),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_smooth(data=subset(ceda4, voteshare_toptwo < 0.5),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_point(data=data.graph.all, 
							 aes(x=mid, y=bin.mean, size=n), pch=1, color="black") + 
		# geom_vline(xintercept=0.5) + 
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013),xlim=c(-ceda_all_cct$bws["h", "left"]-0.01,
		                                              ceda_all_cct$bws["h", "right"]+0.01)) +
		scale_x_continuous(breaks=seq(-0.4, 0.4, .04),
											 limits=c(-ceda_all_cct$bws["h", "left"]-0.01,
											 				 ceda_all_cct$bws["h", "right"]+0.01),
											 labels=seq(0.4, -0.4, -.04)) + 
		theme_minimal() + 
		theme(plot.title = element_text(hjust = 0),
		      legend.position = "bottom", legend.text = element_text(size=4.5),legend.title = element_text(size=4.5),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("CA Local") +
		annotate('text', x = 0.45-0.5, y = .6, 0.05, hjust=0, parse=TRUE, color="black",
						 label = paste('hat(tau)==', round(ceda_all_cct$coef['Conventional', ], 2),
						 							'~(list(', round(ceda_all_cct$ci['Robust', 1], 2),
						 							',', round(ceda_all_cct$ci['Robust', 2], 2), '))'))
)
dev.off()


pdf("figures/rdplot_runanywhere_gender_ceda_binned.pdf", width=3, height=3.5)
(ggplot(ceda4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# female candidates:
		geom_smooth(data=subset(ceda4, voteshare_toptwo > 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_smooth(data=subset(ceda4, voteshare_toptwo < 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_point(data=data.graph.female, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#440154FF") + 
		# male candidates:
		geom_smooth(data=subset(ceda4, voteshare_toptwo > 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_smooth(data=subset(ceda4, voteshare_toptwo < 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_point(data=data.graph.male, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#21908CFF") + 
		# geom_vline(xintercept=0.5) + 
    # scale_color_viridis_d("Candidate gender") + 
    
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013),xlim=c(-ceda_female_cct$bws["h", "left"]-0.0025,
		                                              ceda_female_cct$bws["h", "right"]+0.0025)) +
		scale_x_continuous(breaks=seq(-0.04, 0.04, .02),
											 limits=c(-ceda_female_cct$bws["h", "left"]-0.01,
											 				 ceda_female_cct$bws["h", "right"]+0.01),
											 labels=seq(0.04, -0.04, -.02)) + 
		theme_minimal() + 
		theme(plot.title = element_text(hjust = 0),
		      legend.position = "bottom", legend.text = element_text(size=4.0),legend.title = element_text(size=5.0),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("CA Local") +
		annotate('text', x = -0.057, y = .6, 0.05, hjust=0, parse=TRUE, color="#21908CFF",
						 label = paste('Male:~hat(tau)==', round(ceda_male_cct$coef['Conventional', ], 2),
						 							'~(list(', round(ceda_male_cct$ci['Robust', 1], 2),
						 							',', round(ceda_male_cct$ci['Robust', 2], 2), '))')) + 
		annotate('text', x = -0.057, y = .7, 0.05, hjust=0, parse=TRUE, color="#440154FF",
						 label = paste('Female:~hat(tau)==', round(ceda_female_cct$coef['Conventional', ], 2),
						 							'~(list(', round(ceda_female_cct$ci['Robust', 1], 2),
						 							',', round(ceda_female_cct$ci['Robust', 2], 2), '))')) 
)
dev.off()



# --------------------- #
#### Figures: Mayors ####
# --------------------- #
### Mayors - descriptive
mayors_summ <- mayors4 %>%
  group_by(voteshare>=0.5) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(lose = !win)

plot_descriptive_mayors <- ggplot(mayors_summ, aes(x=lose,y=run_again_anywhere)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.04,label=round(run_again_anywhere*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal()
ggsave(plot_descriptive_mayors,filename = "figures/descriptive_mayors.pdf",width=3,height=3)


plot_descriptive3_mayors <- ggplot(mayors_summ, aes(x=lose,y=run_again_next4yrs)) + 
  geom_bar(stat="identity") + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.04,label=round(run_again_next4yrs*100,0))) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal()
ggsave(plot_descriptive3_mayors,filename = "figures/descriptive3_mayors.pdf",width=3,height=3)

# prep for binned plot:
width <- .01
data.graph <- mayors4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere)))

plot_descriptive_cont_mayors <- ggplot(data.graph.all) + 
  geom_point(aes(x=mid,y=bin.mean,size=n),position=position_jitter(height=0.05),shape=1) + 
  geom_smooth(data=filter(mayors4,voteshare<0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  geom_smooth(data=filter(mayors4,voteshare>=0.5),aes(x=voteshare,y=run_again_anywhere),method="gam") + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal() + 
  theme(legend.position = "bottom")
ggsave(plot_descriptive_cont_mayors,filename = "figures/descriptive_cont_mayors.pdf",width=4,height=4)

## with gender:
width <- .01
data.graph <- mayors4 %>%
  mutate(bin=cut(voteshare, breaks=seq(0.0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere,female) %>%
  group_by(bin,mid,female) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere)))


mayors4 <- mayors4 %>%
  mutate(female2 = recode(female,`1`= "Female",`0`="Male"))
plot_descriptive_cont_nopts_nose_gender_mayors <- ggplot(data.graph.all) + 
  # geom_point(filter(data.graph.all,female==1),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#440154FF",alpha=0.8) + 
  # geom_point(filter(data.graph.all,female==0),aes(x=mid,y=bin.mean,size=n),
  #            position=position_jitter(height=0.05),shape=16,color="#21908CFF",alpha=0.8) + 
  geom_smooth(data=filter(mayors4,voteshare<0.5 & female2=="Female"),aes(x=voteshare,y=run_again_anywhere,col=female2,lty=female2),method="loess",se=F) + 
  geom_smooth(data=filter(mayors4,voteshare>=0.5 & female2=="Female"),aes(x=voteshare,y=run_again_anywhere,col=female2,lty=female2),method="loess",col="#440154FF",se=F) + 
  geom_smooth(data=filter(mayors4,voteshare<0.5 & female2=="Male"),aes(x=voteshare,y=run_again_anywhere,col=female2,lty=female2),method="loess",se=F) + 
  geom_smooth(data=filter(mayors4,voteshare>=0.5 & female2=="Male"),aes(x=voteshare,y=run_again_anywhere,col=female2,lty=female2),method="loess",se=F) + 
  scale_x_continuous("Candidate voteshare",labels=scales::percent_format()) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(-0.005,1.005)) + 
  scale_color_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  scale_linetype_manual("Candidate gender",breaks = c("Female","Male"),values=c(1,2)) + 
  scale_size_continuous("Observations in bin") + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal() + 
  theme(legend.position = c(0.8,0.15))
ggsave(plot_descriptive_cont_nopts_nose_gender_mayors,filename = "figures/descriptive_cont_nopts_nose_gender_mayors.pdf",width=4,height=4)


mayors_summ_gender <- mayors4 %>%
  group_by(voteshare>=0.5,female) %>%
  summarize(run_again_anywhere = mean(run_again_anywhere,na.rm=T),
            run_again_match = mean(run_next_elec_office_match,na.rm=T),
            run_again_next4yrs = mean(run_again_next4yrs,na.rm=T)) %>%
  na.omit() %>%
  rename(win = `voteshare >= 0.5`) %>%
  mutate(female=factor(female,levels = c(0,1),labels=c("Male","Female")),
         lose = !win) %>%
  ungroup()


plot_descriptive_mayors_gender <- ggplot(mayors_summ_gender, aes(x=lose,y=run_again_anywhere,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_anywhere+0.03,label=round(run_again_anywhere*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again anytime",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive_mayors_gender,filename = "figures/descriptive_gender_mayors.pdf",width=6,height=4)



plot_descriptive3_mayors_gender <- ggplot(mayors_summ_gender, aes(x=lose,y=run_again_next4yrs,fill=female)) + 
  geom_bar(stat="identity",position=position_dodge()) + 
  geom_text(aes(x=lose,y=run_again_next4yrs+0.03,label=round(run_again_next4yrs*100,0),group=female),position=position_dodge(width=.9)) + 
  scale_x_discrete("Candidate outcome",breaks=c(TRUE, FALSE),labels=c("Lose","Win")) + 
  scale_y_continuous("% running again in next 4 years",labels=scales::percent_format(),limits=c(0,1)) + 
  # scale_fill_viridis_d("Candidate gender") + 
  scale_fill_manual("Candidate gender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF")) + 
  ggtitle("Nationwide Mayoral") +
  theme_minimal() + 
  theme(legend.position = c(0.75,0.9))
ggsave(plot_descriptive3_mayors_gender,filename = "figures/descriptive3_gender_mayors.pdf",width=6,height=4)



# prep data for binned means:
width <- .005
data.graph <- mayors4 %>%
	mutate(bin=cut(lossmargin, breaks=seq(-0.5,0.5, width))
	)
bins <- data.frame(bin = levels(data.graph$bin),
									 mid = seq(-0.5 + width/2, 0.5 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.female <- data.graph %>% 
	filter(female==1) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin,mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

data.graph.male <- data.graph %>% 
	filter(female==0) %>%
	select(bin, mid, run_again_anywhere) %>%
	group_by(bin, mid) %>%
	summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
						n = sum(!is.na(run_again_anywhere)))

pdf("figures/rdplot_runanywhere_all_mayors_binned.pdf", width=3, height=3.5)
(ggplot(mayors4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# all candidates:
		geom_smooth(data=subset(mayors4, voteshare > 0.5 ),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_smooth(data=subset(mayors4, voteshare < 0.5 ),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="black") + 
		geom_point(data=data.graph.all, 
							 aes(x=mid, y=bin.mean, size=n), pch=1, color="black") + 
		# geom_vline(xintercept=0.5) + 
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013)) +
		scale_x_continuous(breaks=seq(-0.4, 0.4, .04),
											 limits=c(-mayors_all_cct$bws["h", "left"],
											 				 mayors_all_cct$bws["h", "right"]),
											 labels=seq(0.4, -0.4, -.04)) + 
		theme_minimal() + 
		theme(plot.title = element_text(hjust = 0),
		      legend.position = "bottom", legend.text = element_text(size=5.5),legend.title = element_text(size=6.5),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("Nationwide Mayoral") +
		annotate('text', x = 0.45-0.5, y = .85, 0.05, hjust=0, parse=TRUE, color="black",
						 label = paste('hat(tau)==', round(mayors_all_cct$coef['Conventional', ], 2),
						 							'~(list(', round(mayors_all_cct$ci['Robust', 1], 2),
						 							',', round(mayors_all_cct$ci['Robust', 2], 2), '))'))
)
dev.off()


pdf("figures/rdplot_runanywhere_gender_mayors_binned.pdf", width=3, height=3.5)
(ggplot(mayors4) + 
		# geom_point(aes(x = voteshare, y = run_again_anywhere), alpha=.1,pch="l") +  # rugplot slows down plotting
		# female candidates:
		geom_smooth(data=subset(mayors4, voteshare > 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_smooth(data=subset(mayors4, voteshare < 0.5 & female==1),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#440154FF") + 
		geom_point(data=data.graph.female, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#440154FF") + 
		# male candidates:
		geom_smooth(data=subset(mayors4, voteshare > 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_smooth(data=subset(mayors4, voteshare < 0.5 & female==0),
								aes(x = lossmargin, y = run_again_anywhere),
								method = 'lm', formula = y ~ poly(x, 1), size=1, color="#21908CFF",lty=2) + 
		geom_point(data=data.graph.male, 
							 aes(x=mid, y=bin.mean, size=n), pch=16, color="#21908CFF") + 
		# geom_vline(xintercept=0.5) + 
		geom_hline(yintercept=0, lty='dotted') + 
		geom_hline(yintercept=1, lty='dotted') + 
		coord_cartesian(ylim = c(-.013, 1.013),xlim=c(-mayors_female_cct$bws["h", "left"],
		                                             mayors_female_cct$bws["h", "right"]+0.005)) +
		scale_x_continuous(breaks=seq(-0.4, 0.4, .04),
											 limits=c(-mayors_female_cct$bws["h", "left"]-0.01,
											 				 mayors_female_cct$bws["h", "right"]+0.005),
											 labels=seq(0.4, -0.4, -.04)) + 
		theme_minimal() + 
		theme(plot.title = element_text(hjust = 0),
		      legend.position = "bottom", legend.text = element_text(size=5.5),legend.title = element_text(size=6.5),legend.box.margin = unit(c(-1,0,0,-1),"lines")) + 
    labs(size = '# Obs. in\n0.5% Bin',
				 x = "Win Margin",
				 y = "Pr. Running Again Anytime") + 
		ggtitle("Nationwide Mayoral") +
		annotate('text', x = -0.121, y = .9, 0.05, hjust=0, parse=TRUE, color="#21908CFF",
						 label = paste('Male:~hat(tau)==', round(mayors_male_cct$coef['Conventional', ], 2),
						 							'~(list(', round(mayors_male_cct$ci['Robust', 1], 2),
						 							',', round(mayors_male_cct$ci['Robust', 2], 2), '))')) + 
		annotate('text', x = -0.121, y = .1, 0.05, hjust=0, parse=TRUE, color="#440154FF",
						 label = paste('Female:~hat(tau)==', round(mayors_female_cct$coef['Conventional', ], 2),
						 							'~(list(', round(mayors_female_cct$ci['Robust', 1], 2),
						 							',', round(mayors_female_cct$ci['Robust', 2], 2), '))')) 
)
dev.off()


#### Figures with all jurisdictions together ####

## Coefplots:
## Running anytime:

coefs <- rd.export.numeric(list(sler_all_cct,
                                ceda_all_cct,
                                mayors_all_cct
))
coefs$outcome <- c("sler",
                   "ceda",
                   "mayors"
)
coefs_90 <- rd.export.numeric(list(sler_all_cct_90,
                                   ceda_all_cct_90,
                                   mayors_all_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler",
                      "ceda",
                      "mayors"
)

coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "California\nLocal",
                         "Nationwide\nMayoral"
)

pdf("figures/coefplot_anywhere_nogender.pdf",height=3,width=7)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=nrow(coefs):1,ymin=cilo, ymax=cihi), 
                colour="black", width=0, size=0.5) +
  geom_errorbar(aes(x=nrow(coefs):1,ymin=cilo_90, ymax=cihi_90),
                colour="black", width=0, size=1) +
  geom_point(aes(x=nrow(coefs):1,y=coef),size=4, pch=21, fill="black") +
  scale_y_continuous("RD effect of winning vs. losing on probability of\nrunning again anytime + robust CI",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.55,0),
                     labels=seq(-1,0, .1)) +
  scale_x_continuous("Jurisdiction",
                     breaks=c(nrow(coefs):1),
                     limits=c(0.7,nrow(coefs)+0.3),
                     labels=coefs$outcome_pretty,position = "bottom") + 
  coord_flip() + 
  geom_text(aes(x=nrow(coefs):1,y=coef+0.03,label=round(coef,2)),nudge_x=0.06, size=3) +
  theme_bw()
dev.off()

coefs <- rd.export.numeric(list(sler_male_cct,
                                sler_female_cct,
                                ceda_male_cct,
                                ceda_female_cct,
                                mayors_male_cct,
                                mayors_female_cct
))
coefs$outcome <- c("sler-male",
                   "sler-female",
                   "ceda-male",
                   "ceda-female",
                   "mayors-male",
                   "mayors-female"
)

coefs_90 <- rd.export.numeric(list(sler_male_cct_90,
                                   sler_female_cct_90,
                                   ceda_male_cct_90,
                                   ceda_female_cct_90,
                                   mayors_male_cct_90,
                                   mayors_female_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler-male",
                      "sler-female",
                      "ceda-male",
                      "ceda-female",
                      "mayors-male",
                      "mayors-female"
)
coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "Nationwide\nState Legislative",
                         "California\nLocal",
                         "California\nLocal",
                         "Nationwide\nMayoral",
                         "Nationwide\nMayoral"
)
coefs$gender = rep(c("Male","Female"),3)
coefs$order <- c(3,3,2,2,1,1)


pdf("figures/coefplot_anywhere_bygender.pdf",height=4,width=4)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=order,ymin=cilo, ymax=cihi,group=gender,color=gender),position=position_dodge(width = 0.55), 
                width=0, size=1) +
  geom_errorbar(aes(x=order,ymin=cilo_90, ymax=cihi_90,group=gender,color=gender),position=position_dodge(width = 0.55),
                width=0, size=2) +
  geom_point(aes(x=order,y=coef,group=gender,fill=gender,shape=gender,color=gender),position=position_dodge(width = 0.55),size=4) +
  scale_y_continuous("RD effect of losing vs. winning on probability of\nrunning again anytime + robust CI",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.55,0),
                     labels=seq(-1,0, 0.1)) +
  ggtitle("Jurisdiction")+
  scale_x_continuous("",
                     breaks=c(max(coefs$order):1),
                     limits=c(0.7,max(coefs$order)+0.3),
                     labels=unique(coefs$outcome_pretty),position = "bottom") + 
  scale_shape_discrete("Candidate\ngender",solid = T) + 
  # scale_color_viridis_d("Candidate\nGender",aesthetics = c("color","fill")) + 
  scale_color_manual("Candidate\ngender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF"),aesthetics = c("color","fill")) + 
  coord_flip() + 
  geom_text(aes(x=order+0.15,
                y=coef+0.03,
                label=round(coef,2),
                group=gender),
            position=position_dodge(width = 0.55), size=3) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = -0.50),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        legend.position = c(0.25,0.5),
        legend.box.background = element_rect(color = "black"))
dev.off()

## Run in next 4 years
coefs <- rd.export.numeric(list(sler_next4_all_cct,
                                ceda_next4_all_cct,
                                mayors_next4_all_cct
))
coefs$outcome <- c("sler",
                   "ceda",
                   "mayors"
)
coefs_90 <- rd.export.numeric(list(sler_next4_all_cct_90,
                                   ceda_next4_all_cct_90,
                                   mayors_next4_all_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler",
                      "ceda",
                      "mayors"
)

coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "California\nLocal",
                         "Nationwide\nMayoral"
)

pdf("figures/coefplot_next4_nogender.pdf",height=3,width=7)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=nrow(coefs):1,ymin=cilo, ymax=cihi), 
                colour="black", width=0, size=0.5) +
  geom_errorbar(aes(x=nrow(coefs):1,ymin=cilo_90, ymax=cihi_90),
                colour="black", width=0, size=1) +
  geom_point(aes(x=nrow(coefs):1,y=coef),size=4, pch=21, fill="black") +
  scale_y_continuous("RD effect of winning vs. losing on probability of\nrunning again in next 4 years + robust CI",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.55,0),
                     labels=seq(-1,0, .1)) +
  scale_x_continuous("Jurisdiction",
                     breaks=c(nrow(coefs):1),
                     limits=c(0.7,nrow(coefs)+0.3),
                     labels=coefs$outcome_pretty,position = "bottom") + 
  coord_flip() + 
  geom_text(aes(x=nrow(coefs):1,y=coef+0.03,label=round(coef,2)),nudge_x=0.06, size=3) +
  theme_bw()
dev.off()

coefs <- rd.export.numeric(list(sler_next4_male_cct,
                                sler_next4_female_cct,
                                ceda_next4_male_cct,
                                ceda_next4_female_cct,
                                mayors_next4_male_cct,
                                mayors_next4_female_cct
))
coefs$outcome <- c("sler-male",
                   "sler-female",
                   "ceda-male",
                   "ceda-female",
                   "mayors-male",
                   "mayors-female"
)

coefs_90 <- rd.export.numeric(list(sler_next4_male_cct_90,
                                   sler_next4_female_cct_90,
                                   ceda_next4_male_cct_90,
                                   ceda_next4_female_cct_90,
                                   mayors_next4_male_cct_90,
                                   mayors_next4_female_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler-male",
                      "sler-female",
                      "ceda-male",
                      "ceda-female",
                      "mayors-male",
                      "mayors-female"
)
coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "Nationwide\nState Legislative",
                         "California\nLocal",
                         "California\nLocal",
                         "Nationwide\nMayoral",
                         "Nationwide\nMayoral"
)
coefs$gender = rep(c("Male","Female"),3)
coefs$order <- c(3,3,2,2,1,1)


pdf("figures/coefplot_next4_bygender.pdf",height=4,width=4)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=order,ymin=cilo, ymax=cihi,group=gender,color=gender),position=position_dodge(width = 0.55), 
                width=0, size=1) +
  geom_errorbar(aes(x=order,ymin=cilo_90, ymax=cihi_90,group=gender,color=gender),position=position_dodge(width = 0.55),
                width=0, size=2) +
  geom_point(aes(x=order,y=coef,group=gender,fill=gender,shape=gender,color=gender),position=position_dodge(width = 0.55),size=4) +
  scale_y_continuous("RD effect of losing vs. winning on probability of\nrunning again in next 4 years + robust CI",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.55,0),
                     labels=seq(-1,0, 0.1)) +
  ggtitle("Jurisdiction")+
  scale_x_continuous("",
                     breaks=c(max(coefs$order):1),
                     limits=c(0.7,max(coefs$order)+0.3),
                     labels=unique(coefs$outcome_pretty),position = "bottom") + 
  scale_shape_discrete("Candidate\ngender",solid = T) + 
  # scale_color_viridis_d("Candidate\nGender",aesthetics = c("color","fill")) + 
  scale_color_manual("Candidate\ngender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF"),aesthetics = c("color","fill")) + 
  coord_flip() + 
  geom_text(aes(x=order+0.15,
                y=coef+0.03,
                label=round(coef,2),
                group=gender),
            position=position_dodge(width = 0.55), size=3) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = -0.50),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        legend.position = c(0.25,0.5),
        legend.box.background = element_rect(color = "black"))
dev.off()


## First-time candidates only
coefs <- rd.export.numeric(list(sler_firsttime_male_cct,
                                sler_firsttime_female_cct,
                                ceda_firsttime_male_cct,
                                ceda_firsttime_female_cct,
                                mayors_firsttime_male_cct,
                                mayors_firsttime_female_cct
))
coefs$outcome <- c("sler-male",
                   "sler-female",
                   "ceda-male",
                   "ceda-female",
                   "mayors-male",
                   "mayors-female"
)

coefs_90 <- rd.export.numeric(list(sler_firsttime_male_cct_90,
                                   sler_firsttime_female_cct_90,
                                   ceda_firsttime_male_cct_90,
                                   ceda_firsttime_female_cct_90,
                                   mayors_firsttime_male_cct_90,
                                   mayors_firsttime_female_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler-male",
                      "sler-female",
                      "ceda-male",
                      "ceda-female",
                      "mayors-male",
                      "mayors-female"
)
coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "Nationwide\nState Legislative",
                         "California\nLocal",
                         "California\nLocal",
                         "Nationwide\nMayoral",
                         "Nationwide\nMayoral"
)
coefs$gender = rep(c("Male","Female"),3)
coefs$order <- c(3,3,2,2,1,1)


pdf("figures/coefplot_firsttime_bygender.pdf",height=4,width=4)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=order,ymin=cilo, ymax=cihi,group=gender,color=gender),position=position_dodge(width = 0.55), 
                width=0, size=1) +
  geom_errorbar(aes(x=order,ymin=cilo_90, ymax=cihi_90,group=gender,color=gender),position=position_dodge(width = 0.55),
                width=0, size=2) +
  geom_point(aes(x=order,y=coef,group=gender,fill=gender,shape=gender,color=gender),position=position_dodge(width = 0.55),size=4) +
  scale_y_continuous("RD effect of losing vs. winning on probability of\nrunning again anytime\n(first-time candidates only)",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.57,0),
                     labels=seq(-1,0, 0.1)) +
  ggtitle("Jurisdiction")+
  scale_x_continuous("",
                     breaks=c(max(coefs$order):1),
                     limits=c(0.7,max(coefs$order)+0.3),
                     labels=unique(coefs$outcome_pretty),position = "bottom") + 
  scale_shape_discrete("Candidate\ngender",solid = T) + 
  # scale_color_viridis_d("Candidate\nGender",aesthetics = c("color","fill")) + 
  scale_color_manual("Candidate\ngender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF"),aesthetics = c("color","fill")) + 
  coord_flip() + 
  geom_text(aes(x=order+0.15,
                y=coef+0.03,
                label=round(coef,2),
                group=gender),
            position=position_dodge(width = 0.55), size=3) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = -0.50),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        legend.position = c(0.25,0.5),
        legend.box.background = element_rect(color = "black"))
dev.off()


## Run in next 4 years, first-time cands only
coefs <- rd.export.numeric(list(sler_next4_firsttime_male_cct,
                                sler_next4_firsttime_female_cct,
                                ceda_next4_firsttime_male_cct,
                                ceda_next4_firsttime_female_cct,
                                mayors_next4_firsttime_male_cct,
                                mayors_next4_firsttime_female_cct
))
coefs$outcome <- c("sler-male",
                   "sler-female",
                   "ceda-male",
                   "ceda-female",
                   "mayors-male",
                   "mayors-female"
)

coefs_90 <- rd.export.numeric(list(sler_next4_firsttime_male_cct_90,
                                   sler_next4_firsttime_female_cct_90,
                                   ceda_next4_firsttime_male_cct_90,
                                   ceda_next4_firsttime_female_cct_90,
                                   mayors_next4_firsttime_male_cct_90,
                                   mayors_next4_firsttime_female_cct_90
)) %>%
  select(outcome, cilo_90 = cilo, cihi_90 = cihi)
coefs_90$outcome <- c("sler-male",
                      "sler-female",
                      "ceda-male",
                      "ceda-female",
                      "mayors-male",
                      "mayors-female"
)
coefs <- left_join(coefs,coefs_90,by="outcome")

coefs$outcome_pretty = c("Nationwide\nState Legislative",
                         "Nationwide\nState Legislative",
                         "California\nLocal",
                         "California\nLocal",
                         "Nationwide\nMayoral",
                         "Nationwide\nMayoral"
)
coefs$gender = rep(c("Male","Female"),3)
coefs$order <- c(3,3,2,2,1,1)


pdf("figures/coefplot_next4_firsttime_bygender.pdf",height=4,width=4)
ggplot(coefs) +
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(x=order,ymin=cilo, ymax=cihi,group=gender,color=gender),position=position_dodge(width = 0.55), 
                width=0, size=1) +
  geom_errorbar(aes(x=order,ymin=cilo_90, ymax=cihi_90,group=gender,color=gender),position=position_dodge(width = 0.55),
                width=0, size=2) +
  geom_point(aes(x=order,y=coef,group=gender,fill=gender,shape=gender,color=gender),position=position_dodge(width = 0.55),size=4) +
  scale_y_continuous("RD effect of losing vs. winning on probability of\nrunning again in next 4 years\n(first-time candidates only)",
                     breaks=seq(-1,0, 0.1),
                     limits=c(-0.57,0),
                     labels=seq(-1,0, 0.1)) +
  ggtitle("Jurisdiction")+
  scale_x_continuous("",
                     breaks=c(max(coefs$order):1),
                     limits=c(0.7,max(coefs$order)+0.3),
                     labels=unique(coefs$outcome_pretty),position = "bottom") + 
  scale_shape_discrete("Candidate\ngender",solid = T) + 
  # scale_color_viridis_d("Candidate\nGender",aesthetics = c("color","fill")) + 
  scale_color_manual("Candidate\ngender",breaks=c("Female","Male"),values = c("#440154FF","#21908CFF"),aesthetics = c("color","fill")) + 
  coord_flip() + 
  geom_text(aes(x=order+0.15,
                y=coef+0.03,
                label=round(coef,2),
                group=gender),
            position=position_dodge(width = 0.55), size=3) +
  theme_bw() + 
  theme(plot.title = element_text(hjust = -0.50),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        legend.position = c(0.25,0.5),
        legend.box.background = element_rect(color = "black"))
dev.off()



# ------------------------------------------- #
#### Tables with all jurisdictions together ####
# ------------------------------------------- #
## RDD results
rdd_tab <- rd.export4(list(sler_male_cct,
                           sler_female_cct,
                           ceda_male_cct,
                           ceda_female_cct,
                           mayors_male_cct,
                           mayors_female_cct)) %>% as.data.frame()
rdd_tab$Jurisdiction <- c("Nationwide State Legislative","","","","CA Local","","","","Nationwide Mayoral","","","")
rdd_tab$Subset <- c("Male","","Female","","Male","","Female","","Male","","Female","")
rdd_tab <- select(rdd_tab,
                  Jurisdiction,Subset,Coef,`p-value`,BW,Obs)

print(xtable(rdd_tab),file = "tables/rdd_bygender.tex",floating = F,include.rownames = F)


rdd_tab <- rd.export4(list(sler_next4_male_cct,
                           sler_next4_female_cct,
                           ceda_next4_male_cct,
                           ceda_next4_female_cct,
                           mayors_next4_male_cct,
                           mayors_next4_female_cct)) %>% as.data.frame()
rdd_tab$Jurisdiction <- c("Nationwide State Legislative","","","","CA Local","","","","Nationwide Mayoral","","","")
rdd_tab$Subset <- c("Male","","Female","","Male","","Female","","Male","","Female","")
rdd_tab <- select(rdd_tab,
                  Jurisdiction,Subset,Coef,`p-value`,BW,Obs)

print(xtable(rdd_tab),file = "tables/rdd_next4_bygender.tex",floating = F,include.rownames = F)


rdd_tab <- rd.export4(list(sler_firsttime_male_cct,
                           sler_firsttime_female_cct,
                           ceda_firsttime_male_cct,
                           ceda_firsttime_female_cct,
                           mayors_firsttime_male_cct,
                           mayors_firsttime_female_cct)) %>% as.data.frame()
rdd_tab$Jurisdiction <- c("Nationwide State Legislative","","","","CA Local","","","","Nationwide Mayoral","","","")
rdd_tab$Subset <- c("Male","","Female","","Male","","Female","","Male","","Female","")
rdd_tab <- select(rdd_tab,
                  Jurisdiction,Subset,Coef,`p-value`,BW,Obs)

print(xtable(rdd_tab),file = "tables/rdd_firsttime_bygender.tex",floating = F,include.rownames = F)


rdd_tab <- rd.export4(list(sler_next4_firsttime_male_cct,
                           sler_next4_firsttime_female_cct,
                           ceda_next4_firsttime_male_cct,
                           ceda_next4_firsttime_female_cct,
                           mayors_next4_firsttime_male_cct,
                           mayors_next4_firsttime_female_cct)) %>% as.data.frame()
rdd_tab$Jurisdiction <- c("Nationwide State Legislative","","","","CA Local","","","","Nationwide Mayoral","","","")
rdd_tab$Subset <- c("Male","","Female","","Male","","Female","","Male","","Female","")
rdd_tab <- select(rdd_tab,
                  Jurisdiction,Subset,Coef,`p-value`,BW,Obs)

print(xtable(rdd_tab),file = "tables/rdd_next4_firsttime_bygender.tex",floating = F,include.rownames = F)




# ------------------- #
#### McCrary tests ####
# ------------------- #

#### SLERs:
## sample one candidate from every election:
nested.sler <- sler4 %>%
  filter(!is.na(voteshare_toptwo)) %>% # only cands in top two with valid voteshare
  group_by(elec_index) %>% # group by election
  nest() %>% # nest within groups
  ungroup()
dim(nested.sler) # 141,116 elections

set.seed(02139)
sample_data_sler <- nested.sler %>%
  mutate(samp = map(data,sample_n,size=1)) %>%
  select(-data) %>%
  unnest(samp)
dim(sample_data_sler) # 141116 candidate-elections



# formal McCrary test:
width <- .005
data.graph <- sample_data_sler %>%
  mutate(bin=cut(voteshare_toptwo, breaks=seq(0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere))) %>%
  mutate(mid_adj = mid - 0.5)

mc.rd.sler <- lm(n ~ mid_adj * (mid_adj>=0), data=data.graph.all[which(data.graph.all$mid>=(0.5-sler_all_cct$bws["h", "left"]) & data.graph.all$mid<=(0.5+sler_all_cct$bws["h", "right"])),])
summary(mc.rd.sler)
stargazer(mc.rd.sler,out = "tables/mccrary.tex",
          dep.var.labels = "Number of observations in bin",
          column.labels = c("Nationwide State Legislative"),
          float = F)

### Density histogram with sampled data:
mccrary_sler <- ggplot(sample_data_sler,aes(x=voteshare_toptwo)) +
  geom_histogram(alpha=0.8,bins = 16) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.02),
                     limits=c(0.44,0.56),labels=scales::percent_format(accuracy=1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("Nationwide State Legislative")

ggsave(mccrary_sler,filename = "figures/mccrary_sample_sler.pdf",height=3,width=4)

mccrary_full_sler <- ggplot(sample_data_sler,aes(x=voteshare_toptwo)) +
  geom_histogram(alpha=0.8,bins = 32) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.2),
                     limits=c(0,1),labels=scales::percent_format(accuracy=1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("Nationwide State Legislative")

ggsave(mccrary_full_sler,filename = "figures/mccrary_sample_full_sler.pdf",height=3,width=4)


#### CEDA:
## sample one candidate from every election:
nested.ceda <- ceda4 %>%
  filter(!is.na(voteshare_toptwo)) %>% # only cands in top two with valid voteshare
  group_by(raceid_new) %>% # group by election
  nest() %>% # nest within groups
  ungroup()
dim(nested.ceda) # 21,875 elections

set.seed(02139)
sample_data_ceda <- nested.ceda %>%
  mutate(samp = map(data,sample_n,size=1)) %>%
  select(-data) %>%
  unnest(samp)
dim(sample_data_ceda) # 21,875 candidate-elections

# formal McCrary test:
width <- .005
data.graph <- sample_data_ceda %>%
  mutate(bin=cut(voteshare_toptwo, breaks=seq(0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere))) %>%
  mutate(mid_adj = mid - 0.5)

mc.rd.ceda <- lm(n ~ mid_adj * (mid_adj>=0), data=data.graph.all[which(data.graph.all$mid>=(0.5-ceda_all_cct$bws["h", "left"]) & data.graph.all$mid<=(0.5+ceda_all_cct$bws["h", "right"])),])
summary(mc.rd.ceda)
stargazer(mc.rd.sler,mc.rd.ceda,
          out = "tables/mccrary.tex",
          dep.var.labels = "Number of observations in bin",
          column.labels = c("Nationwide State Legislative","CA Local"),
          float = F)

### Density histogram with sampled data:
mccrary_ceda <- ggplot(sample_data_ceda,aes(x=voteshare_toptwo)) +
  geom_histogram(alpha=0.8, bins=16) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.02),
                     limits=c(0.44,0.56),labels=scales::percent_format(accuracy=1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("CA Local")

ggsave(mccrary_ceda,filename = "figures/mccrary_sample_ceda.pdf",height=3,width=4)

mccrary_full_ceda <- ggplot(sample_data_ceda,aes(x=voteshare_toptwo)) +
  geom_histogram(alpha=0.8, bins=32) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.2),
                     limits=c(0,1),labels=scales::percent_format(accuracy=1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("CA Local")

ggsave(mccrary_full_ceda,filename = "figures/mccrary_sample_full_ceda.pdf",height=3,width=4)


#### Mayors:
## sample one candidate from every election:
nested.mayors <- mayors4 %>%
  filter(!is.na(voteshare)) %>% # only cands in top two with valid voteshare
  group_by(elec_index) %>% # group by election
  nest() %>% # nest within groups
  ungroup()
dim(nested.mayors) # 7714 elections

set.seed(02139)
sample_data_mayors <- nested.mayors %>%
  mutate(samp = map(data,sample_n,size=1)) %>%
  select(-data) %>%
  unnest(samp)
dim(sample_data_mayors) # 21,875 candidate-elections

# formal McCrary test:
width <- .005
data.graph <- sample_data_mayors %>%
  mutate(bin=cut(voteshare, breaks=seq(0,1, width))
  )
bins <- data.frame(bin = levels(data.graph$bin),
                   mid = seq(0 + width/2, 1 - width/2, width)
)
data.graph <- left_join(data.graph,bins,by="bin")

data.graph.all <- data.graph %>% 
  select(bin, mid, run_again_anywhere) %>%
  group_by(bin,mid) %>%
  summarise(bin.mean = mean(run_again_anywhere,na.rm=T),
            n = sum(!is.na(run_again_anywhere))) %>%
  mutate(mid_adj = mid - 0.5)

mc.rd.mayors <- lm(n ~ mid_adj * (mid_adj>=0), data=data.graph.all[which(data.graph.all$mid>=(0.5-mayors_all_cct$bws["h", "left"]) & data.graph.all$mid<=(0.5+mayors_all_cct$bws["h", "right"])),])
summary(mc.rd.mayors)

stargazer(mc.rd.sler, mc.rd.ceda, mc.rd.mayors,
          out = "tables/mccrary.tex",
          dep.var.labels = "Number of observations in bin",
          column.labels = c("Nationwide State Legislative","CA Local","Nationwide Mayoral"),
          omit.stat = c("ser","f","adj.rsq"),
          covariate.labels = c("Voteshare bin","Voteshare $\\ge 0.5$","Voteshare bin $\\times$ Voteshare $\\ge$ 0.5"),
          float = F)

### Density histogram with sampled data:
mccrary_mayors <- ggplot(sample_data_mayors,aes(x=voteshare)) +
  geom_histogram(alpha=0.8, bins=16) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.02),
                     limits=c(0.44,0.56),labels=scales::percent_format(accuracy = 1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("Nationwide Mayoral")

ggsave(mccrary_mayors,filename = "figures/mccrary_sample_mayors.pdf",height=3,width=4)

mccrary_full_mayors <- ggplot(sample_data_mayors,aes(x=voteshare)) +
  geom_histogram(alpha=0.8, bins=32) +
  geom_vline(col="red",lty=2,xintercept = 0.5) +
  labs(x="Voteshare",y="Observations") +
  scale_x_continuous(breaks=seq(0,1,0.2),
                     limits=c(0,1),labels=scales::percent_format(accuracy = 1)) +
  theme(text = element_text(size=20)) + 
  theme_minimal() + 
  ggtitle("Nationwide Mayoral")

ggsave(mccrary_full_mayors,filename = "figures/mccrary_sample_full_mayors.pdf",height=3,width=4)


## Using Cattaneo Jansson and Ma (2019) nonparametric method:
library(rddensity)
cjm_sler <- with(sample_data_sler, rddensity(X=lossmargin,c=0,massPoints = F))
summary(cjm_sler)

cjm_ceda <- with(sample_data_ceda, rddensity(X=lossmargin,c=0,massPoints = F))
summary(cjm_ceda)

cjm_mayors <- with(sample_data_mayors, rddensity(X=lossmargin,c=0,massPoints = F))
summary(cjm_mayors)

names(cjm_sler)
cjm_sler$hat$diff # density difference on either side of cutoff
cjm_sler$test$t_jk # t-stat with SE based on jackknife
cjm_sler$test$p_jk # p-value for desntiy test
cjm_sler$N$eff_left + cjm_sler$N$eff_right # effective N
cjm_sler$h$left
cjm_sler$h$right

density_tab <- data.frame("Jurisdiction" = c("State Legislative","CA Local","Nationwide Mayoral"),
                          "t-statistic" = c(cjm_sler$test$t_jk,cjm_ceda$test$t_jk,cjm_mayors$test$t_jk),
                          "p-value" = c(cjm_sler$test$p_jk,cjm_ceda$test$p_jk,cjm_mayors$test$p_jk),
                          "Effective N" = c(cjm_sler$N$eff_left + cjm_sler$N$eff_right,
                                            cjm_ceda$N$eff_left + cjm_ceda$N$eff_right,
                                            cjm_mayors$N$eff_left + cjm_mayors$N$eff_right)#,
                          # "Bandwidth" = c(cjm_sler$h$left + cjm_sler$h$right,
                          #                 cjm_ceda$h$left + cjm_ceda$h$right,
                          #                 cjm_mayors$h$left + cjm_mayors$h$right)
                          )

print(xtable(density_tab),file = "tables/rddensity_cjm.tex",floating = F,include.rownames = F)


# ----------------------------------- #
#### Null hypothesis of difference ####
# ----------------------------------- #

## For RDD Density test using Hartman (2020) method
devtools::source_url("https://github.com/ekhartman/rdd_equivalence/blob/master/RDD_equivalence_functions.R?raw=TRUE")
fstar_sler <- cjm_sler$hat$left/cjm_sler$hat$right # f-star-hat
hartman_sler <- rdd.tost.ratio(estL = cjm_sler$hat$left, estR = cjm_sler$hat$right,
                               seL = cjm_sler$sd_jk$left, seR = cjm_sler$sd_jk$right,
                               eps = 1.5,
                               alpha = 0.05)
fstar_ceda <- cjm_ceda$hat$left/cjm_ceda$hat$right # f-star-hat
hartman_ceda <- rdd.tost.ratio(estL = cjm_ceda$hat$left, estR = cjm_ceda$hat$right,
                               seL = cjm_ceda$sd_jk$left, seR = cjm_ceda$sd_jk$right,
                               eps = 1.5,
                               alpha = 0.05)
fstar_mayors <- cjm_mayors$hat$left/cjm_mayors$hat$right # f-star-hat
hartman_mayors <- rdd.tost.ratio(estL = cjm_mayors$hat$left, estR = cjm_mayors$hat$right,
                               seL = cjm_mayors$sd_jk$left, seR = cjm_mayors$sd_jk$right,
                               eps = 1.5,
                               alpha = 0.05)


hartman_tab <- data.frame("Jurisdiction" = c("Nationwide State Legislative","CA Local","Nationwide Mayoral"),
                          "Observed Ratio" = c(fstar_sler,fstar_ceda,fstar_mayors),
                          "Equivalence Confidence Interval" = c(paste0("(",round(1/hartman_sler$inverted,2), ", ",round(hartman_sler$inverted,2),")"),
                                                                paste0("(",round(1/hartman_ceda$inverted,2), ", ",round(hartman_ceda$inverted,2),")"),
                                                                paste0("(",round(1/hartman_mayors$inverted,2), ", ",round(hartman_mayors$inverted,2),")")),
                          "p-value" = c(hartman_sler$p,hartman_ceda$p,hartman_mayors$p)
                          )
print(xtable(hartman_tab),file = "tables/density_hartman.tex",floating = F,include.rownames = F)



